!Function ghost_energy(label)
!use type
!use var
!implicit none
!integer :: i,j,k,label
!real(dp) :: LdG = 0.d0, Elastic = 0.d0
!real(dp) :: ghost_energy, acu 
!
!LdG     = LdG_Energy(q(:,label))
!Elastic = Elastic_Energy(label)
!ghost_energy = Elastic + LdG  
!
!return
!contains

Function LdG_Energy(q_tem)
use type
use var
implicit none
real(dp), dimension(6) :: q_tem
real(dp) :: tr1, tr2, tr3
real(dp) :: aux1, a1, a4
real(dp) ::  LdG_Energy

tr1 =  2*q_tem(1)**2 + 2*q_tem(2)**2 + 2*q_tem(3)**2 + 2*q_tem(4)**2 + 2*q_tem(5)**2 + 2*q_tem(1)*q_tem(4)

tr2 = tr1 * tr1

aux1 = 3.d0/2.d0
a1 = -q_tem(1) - q_tem(4)
a4 = q_tem(1) - q_tem(4)

tr3 = 0.75d0*a1**3 - 3.d0*a1*q_tem(2)**2 + aux1*a1*q_tem(3)**2 + aux1*q_tem(3)**2*a4 &
       - 0.5d0*aux1*a1*a4**2 + 6.d0*q_tem(2)*q_tem(3)*q_tem(5) + aux1*a1*q_tem(5)**2 - aux1*a4*q_tem(5)**2
LdG_Energy = 0.5d0*( 1.d0 - ULdG/3.d0 )*tr1 - ULdG/3.d0*tr3 + ULdG/4.d0*tr2
return
end function LdG_Energy
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Function  Elastic_Energy(label_c)
use type
use var
implicit none
integer :: ii
integer :: x, y, left, right, label_l, label_r, label_c
real(dp), dimension(36) :: grad
real(dp) :: Elastic_Energy

Do x = 1 , 3
   ii = (x-1)*2 + 1
   label_l = neighbors(ii,label_c)
   label_r = neighbors(ii+1,label_c)
   grad(1+(x-1)*12:x*12-6) = (q(:,label_r)-q(:,label_c))*df(x) 
   grad(7+(x-1)*12:x*12)   = (q(:,label_c)-q(:,label_l))*df(x)
   grad(1+(x-1)*12:x*12-6) = grad(1+(x-1)*12:x*12-6)*grad(1+(x-1)*12:x*12-6)*fac
   grad(7+(x-1)*12:x*12) = grad(7+(x-1)*12:x*12)*grad(7+(x-1)*12:x*12)*fac
enddo

Elastic_Energy = 0.5d0*sum(grad)
return
end function Elastic_Energy

!End function  ghost_energy


